home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 007 / powern.arc / POWERN.INC
Text File  |  1985-07-22  |  2KB  |  64 lines
  1. {       POWERN.PLB #1.04 85-06-22 FAST INTEGER POWERS OF REAL NUMBERS
  2.  
  3.       V01 L04  pretty-print update on 85-06-22 with expanded remarks by DEH
  4.           L03  update on 85-06-18 by DEH to expand on an unpublished idea
  5.                of David Gries.
  6.           L02  update on 85-06-11 by DEH for draft DDJ pretty printing.
  7.           L01  update on 85-06-09 by DEH to smooth out some rough edges.
  8.           L00  created on 85-06-05 by Dennis E. Hamilton in order to answer
  9.                a number of electronic requests for help with Pascal powers.}
  10.  
  11. function
  12.  
  13.    PowerN(x: real;   {initial value: x0}
  14.           n: integer {required exponent} )
  15.  
  16.        :real {value of x0^n or x0**n};
  17.  
  18.    var    i: integer {non-negative power of x which remains to be computed};
  19.           r: real {running intermediate-result value};
  20.  
  21.    procedure
  22.  
  23.       MakeOdd;
  24.  
  25.       begin {V01 L03 Knuth transformation re-ordered to exploit
  26.              initial i > 0 and even, thanks to the Gries adjustment.}
  27.       repeat {reduction of i while preserving x0^abs(n) = r*x^i}
  28.          x := sqr(x);
  29.          i := i shr 1;
  30.               {or use i := i div 2 when shifting is not supported}
  31.          until odd(i);
  32.       end; {transformation of x^(2*k) to sqr(x)^k until odd k reached}
  33.  
  34.    BEGIN {PowerN}
  35.  
  36.    if n = 0
  37.    then begin
  38.         if abs(x) = 0.0
  39.                     {V01 L01: Insulating against Turbo Pascal and
  40.                      IEEE -0.0 possibilities}
  41.         then PowerN := x/x
  42.                   {In the indeterminate 0^0 case, defer to the system's
  43.                    treatment of 0/0 for appropriate continuation.}
  44.         else PowerN := 1.0;
  45.         end
  46.    else begin {n <> 0 now guaranteed}
  47.         i := abs(n);
  48.         if not odd(i) then MakeOdd;
  49.           {V01 L03 odd(i) now guaranteed}
  50.         r := x; i := pred(i);
  51.           {V01 L03 even(i) now guaranteed for the Gries invariant}
  52.         while i <> 0
  53.            do begin {not odd(i) and x0^abs(n) = r*x^i}
  54.               MakeOdd;
  55.               r := r*x; i := pred(i);
  56.               end {Gries-invariant version of the Knuth transformation
  57.                    as re-expressed in Dijkstra-Jensen-Wirth form};
  58.         if n < 0
  59.         then PowerN := 1.0/r
  60.         else PowerN := r;
  61.         end;
  62.  
  63.    END {PowerN};
  64.